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Successful  satellite  operations  often  depend  on  the  ability  to  precisely  determine  the 
satellite’s  position  throughout  it’s  operational  lifetime.  Often,  simplified  models  are  used  to 
reduce  the  computational  complexity  when  producing  an  orbit  determination  solution.  One 
example  is  the  use  of  lower  order  gravitational  terms  when  computing  the  state  transition 
matrix.  The  US  Naval  Research  Laboratory’s  Orbit  Covariance  Estimation  and  Analysis 
software  is  used  to  characterize  the  affects  of  using  the  full  order  geopotential  terms  when 
calculating  the  state  transition  matrix.  Orbit  solutions  are  calculated  for  several  geodetic 
satellite  using  satellite  laser  ranging  data.  The  RMS  of  the  error  residual  and  number  of 
iterations  needed  for  convergence  are  used  as  metric  to  evaluate  the  difference  of  using  the 
full  order  geopotential  to  calculate  the  state  transition  matrix.  It  is  found  that  typically  the 
RMS  error  is  the  same  when  the  higher  order  gravitational  terms  are  used;  however,  the 
number  of  iterations  needed  for  convergence  are  often  lower.  Additionally,  the  difference 
in  the  orbit  solution  between  both  methods  is  examined.  Surprisingly,  these  orbit  solutions 
may  differ  by  up  to  several  meters.  This  indicates  that  the  higher  order  gravitational  terms 
should  always  be  considered  when  calculating  the  state  transition  matrix.  Because  fewer 
iterations  are  often  needed,  this  can  result  in  an  overall  decrease  in  computational  time 
required. 


Nomenclature 


x  Satellite  state  vector 
3?"  n-dimensional  real  space 
<1?  Satellite  state  transition  matrix 
State  covariance  matrix 
A  time  varying  dynamics  matrix 
V  the  geopotential 
r  the  Hessian  of  the  geopotential 

J2  the  Earth’s  second  zonal  spherical  harmonic  gravity  term 


I.  Introduction 

Successful  satellite  operations  are  often  dependent  on  the  ability  to  precisely  determine  the  satellite’s 
position  throughout  the  operational  lifetime.  Often,  simplified  models  are  used  to  reduce  the  computational 
complexity  without  sacrificing  fidelity  of  the  precision  orbit  determination  solution.  One  such  simplification 
is  the  use  of  only  the  J2  order  gravitational  terms  when  computing  the  state  transition  matrix.  Using  the 
Orbit  Covariance  Estimation  and  ANalysis  (OCEAN)  precision  orbit  determination  software,  the  affects  of 
using  the  full  order  geopotential  terms  when  calculating  the  state  transition  matrix  is  evaluated. 

Precision  orbit  solutions  are  calculated  for  the  geodetic  satellites  LAGEOS-1,  LAGEOS-2,  STELLA 
and  STARLETTE  using  satellite  laser  ranging  data.  The  orbit  of  these  satellites  is  well  suited  for  orbit 
determination  studies  since  they  have  simple  shapes  and  are  passive  in  nature.  By  design,  their  orbits  are 
dominated  by  gravitational  forces.  So,  error  in  an  orbit  solution  can  be  attributed  to  mis-modeling  forces 
rather  than  to  mis-modeling  satellite  attitude  or  maneuvers. 
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The  RMS  of  the  residual  error  for  a  30  day  fit  span  as  well  as  the  number  of  iterations  required  for 
convergence  will  be  used  as  metrics.  These  metrics  will  allow  for  the  characterization  of  the  convergence 
behavior  as  well  as  the  general  quality  of  the  orbit  solution.  Further,  the  resulting  solution  ephemerides  are 
compared.  This  comparison  allows  for  the  characterization  of  the  differences  in  the  orbit  solution. 

II.  OCEAN  Orbit  Determination  Software 

The  Orbit  Covariance  Estimation  and  Analysis  (OCEAN)  software  was  developed  by  the  Naval  Center 
for  Space  Technology  at  the  Naval  Research  Laboratory  (NRL)  to  not  only  provide  a  tool  to  perform  astro- 
dynamics  related  research,  but  to  also  perform  precision  orbit  determination  to  support  satellite  operations. 

Early  history  of  OCEAN  is  given  in  Reference  1.  References  2,  3  and  4  discuss  further  developments 
of  OCEAN.  More  recently,  OCEAN  has  been  used  to  calculate  orbits  to  support  operations  for  the  NRL 
UPPERSTAGE  and  TACSAT-4  satellites. 

OCEAN  is  a  highly  configurable,  database  driven  software  tool  that  may  be  used  to  calculate  precision 
orbit  solutions  for  a  range  of  satellite  missions.  OCEAN  allows  users  to  simulate  data,  propagate  a  space¬ 
craft  state,  or  solve  for  an  orbit  using  a  Kalman  Filter-Smoother  (KFS)  or  Weighted  Least  Squares  Orbit 
Determination  (WLS-OD)  process. 

III.  Orbit  Determination  Methodology  and  Implementation 

Precision  Orbit  Determination  (POD)  is  often  accomplished  by  using  an  Extended  Kalman  Filter  (EKF) 
or  a  batch  Weighted  Least  Squares  (WLS)  process.5  In  either  process,  the  satellite  dynamics  are  linearized 
around  an  a  priori  solution.  For  both  the  EKF  and  WLS  process,  the  covariance  matrix  is  propagated 
through  time  using  the  state  transition  matrix.  Typically,  the  state  transition  matrix  is  calculated  by 
numerically  integrating  the  dynamics  matrix6.' 

Solving  for  the  state  transition  matrix  is  equivalent  to  solving  initial-value  problems  beginning  with: 

Xj(f)  =  Axj(t),  t  G  3?,  x,(t0)  =  e.i  £  S"  (1) 

Due  to  this  system  being  time- variant,  the  state  transition  matrix  is  defined  as  the  following  for  all  t,  to 
G  3?  : 

*(M.)  =  =  £  A‘V°)t  (2) 

k= 0 

Where  to  and  t  are  the  start  and  end  time  respectively,  $  is  the  state  transition  matrix,  and  A  is  the 
time  varying  dynamics  matrix  in  the  general,  linearized  state  space  model.  For  the  derivation  of  this  state 
transition  matrix  and  related  results,  see  Reference  8. 

The  following  equations  for  the  state  transition  matrix  are  true. 

3>  =  A$,  *(t„)=I  (3) 

This  allows  the  state  transition  matrix  to  be  numerically  integrated  using  the  dynamics  matrix. 

In  the  case  of  the  EKF,  the  covariance  is  used  to  determine  the  Kalman  gain,  which  is  used  to  calculate 
the  state  update.  In  the  case  of  the  WLS  process,  the  covariance  at  each  measurement  is  propagated  to 
the  estimation  epoch  where  it  is  then  used  in  the  normal  equations  to  calculate  the  state  estimate.  The 
covariance  is  propagated  through  time  using  the  following  equation. 

P  fc  =  $(4,4-i)Pfc-i$T(4,4-i)  (4) 

Here,  is  the  covariance  at  time  step  k  and  $>(tk,tk- i)  is  the  state  transition  matrix  from  time  step 
k  —  1  to  time  step  k. 

By  using  higher  fidelity  modeling  the  dynamics  matrix,  a  more  precise  state  transition  matrix  is  calculated 
and  therefore  a  more  precise  covariance  is  calculated  throughout  time.  Reference  9  gives  more  detail  on  the 
mathematics  of  an  EKF  and  WLS  process  as  applied  to  orbit  determination. 

To  demonstrate  how  the  dynamics  matrix  relates  to  the  covariance  matrix,  first  consider  the  equations 
of  motion  for  a  satellite  acted  on  by  only  gravitational  forces. 
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r  =  g  =  W 


(5) 


where  V  is  the  geopotential.10 

For  a  state  x  =  (x,  y,  z,  x,  y,  i),  the  dynamics  matrix  A  is  dehned  in  Reference  11  and  summarized  in 
Reference  10. 


where 


x  =  Ax 

/  0  0  0  1  0  o\ 

0  0  0  0  1  0 
A_  0  0  0001 

rx,x  r y>x  rz>x  o  o  o 
r x,y  Tyty  tZ}V  ooo 
\Tx,z  r VtZ  rZtZ  ooo/ 

d2V  ^  d2V  ^  d2V 


(6) 

(7) 


(8) 


r  represents  the  components  of  the  gravity  gradient  matrix,  or  Hessian  of  the  geopotential,  in  the  rect¬ 
angular  inertial  coordinates.10 

As  a  simplification,  the  Hessian  of  the  geopotential  is  often  calculated  by  considering  lower  order  gravity 
terms  (just  J-i)  as  demonstrated  in  Reference  12.  However,  Reference  10  derives  the  expressions  for  the 
full  order  Hessian  of  the  geopotential  in  terms  of  the  solution  to  Laplace’s  equation.  This  allows  the  state 
transition  matrix  to  be  integrated  using  arbitrary  order  geopotential  model.  Both  the  method  from  Reference 
12  and  Reference  10  are  implemented  within  OCEAN  and  will  be  used  to  demonstrate  impact  of  these 
methods  on  the  WLS  orbit  determination  process. 

Once  data  is  processed  computationally,  it  will  be  compared  to  SLR  results  provided  by  the  International 
Laser  Ranging  Service  (ILRS).  ILRS  provides  this  data  for  several  retro- reflector  fitted  spacecraft  in  order 
to  support  geodetic  research  activities.  Some  of  the  supported  scientific  and  operational  data  products 
include  Earth  Orientation  Parameters,  Fundamental  Physical  Constants,  Monitoring  Earth  Rotation  and 
Polar  Motion,  and  Time- Varying  Geocenter  Coordinates.  For  a  full  list  of  capabilities  refer  to  Reference  13. 

Orbit  solutions  are  calculated  using  the  OCEAN  WLS-OD  process  for  30  days  of  data.  This  long  duration 
orbit  fit  demonstrates  how  well  the  long  term  satellite  dynamics  can  be  modeled.  The  RMS  residual  error 
is  used  as  a  metric  to  demonstrate  how  well  the  data  are  fit.  Further,  the  number  of  iterations  needed  for 
convergence  are  tracked.  This  allows  the  effects  on  the  solution  convergence  for  each  method  to  be  tracked. 
Last,  orbit  ephemeris  solutions  are  compared  to  determine  the  effects  of  both  methods  on  the  solution. 

The  OCEAN  WLS-OD  capability  employs  extensive  spacecraft,  measurement,  and  force  modeling  to 
estimate  the  desired  spacecraft  state  and  parameters.  Force  models  include  solid  Earth  tides;  pole  tides;  lunar 
and  solar  third  body  gravitational  effects;  indirect  lunar  oblateness;  general  relativistic  effects;  atmospheric 
variability;  drag;  and  solar  radiation  pressure.  OCEAN  WLS-OD  accounts  for  various  time  systems,  including 
TAI,  UTC,  and  GPS  time.  The  WLS-OD  functionality  also  uses  Earth  Orientation  Parameter  (EOP)  models 
from  the  International  Earth  Rotation  and  Reference  Systems  Service  (IERS),  which  account  for  precession 
and  nutation,  Earth  rotation,  and  polar  motion.  OCEAN  implements  the  majority  of  the  conventions 
recommended  in  IERS  Conventions  2010. 14  For  this  application,  a  ninth  order  multi-step  predictor-corrector 
algorithm  is  used  to  perform  the  integration  of  the  state  variables  and  state  transition  matrix.  OCEAN  will 
be  used  with  SLR  data  to  calculate  precision  orbits  for  several  geodetic  satellites,  including,  LAGEOS-1, 
LAGEOS-2,  STARLETTE  and  STELLA.  These  satellites  are  well  suited  for  orbit  determination  analysis  as 
their  orbit  was  designed  to  be  sensitive  only  to  the  well  characterized  gravity  forces  and  minimally  sensitive 
to  atmospheric  drag  and  radiation  pressure  forces.  Details  for  each  satellite  are  as  follows. 


A.  LAGEOS-1 

The  LAGEOS-1  satellite  was  launched  in  1976  by  the  National  Aeronautics  and  Space  Administration 
(NASA)  as  part  of  the  Earth  and  ocean  dynamics  application  program.  It  was  designed  to  provide  a  long- 
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lasting  laser  target  in  a  well  defined  orbit.  This  satellite  enabled  researchers  to  study  a  range  of  geophysical 
phenomena  with  improved  accuracy,  including  the  Earth’s  geopotential. 

The  satellite’s  low  ballistic  coefficient  combined  with  its  spherical  shape  minimize  the  orbital  uncertainty 
due  to  drag  and  solar  radiation  forces. 

The  nominal  orbital  elements  are  given  in  Table  1. 


Table  1.  Nominal  Orbital  Elements  for  LAGEOS-1 


Element 

Nominal  Value 

Semi-major  Axis 

Eccentricity 

Inclination 

12,  240  km 

0.0045 

109.84° 

B.  LAGEOS-2 

The  LAGEOS-2  satellite  was  launched  in  1992  by  NASA  and  the  Agenzia  Spaziale  Italiana  (ASI).  The 
LAGEOS-2  satellite  is  an  almost  identical  design  to  the  LAGEOS-1  satellite;  however,  its  orbit  was  selected 
to  provide  increased  coverage  over  seismically  active  areas  such  as  the  Mediterranean  and  California. 

The  nominal  orbital  elements  for  LAGEOS-2  are  given  in  Table  2. 


Table  2.  Nominal  Orbital  Elements  for  LAGEOS-2 


Element 

Nominal  Value 

Semi-major  Axis 

Eccentricity 

Inclination 

1 1 ,  998  km 

0.0135 

52.64° 

C.  STARLETTE 

The  STARLETTE  satellite  was  launched  in  1975  by  Centre  Nationale  d’Etudes  Saptiales  (CNES).  The 
spacecraft  was  designed  to  improve  the  geopotential  model  and  to  study  solid  Earth  tides,  ocean  tides,  and 
polar  motion.15  STARLETTE  was  also  the  first  spacecraft  to  be  entirely  covered  by  laser  corner  reflectors 
which  allow  passive  SLR  observation  capabilities.16  The  STARLETTE  orbit  is  highly  sensitive  to  zonal 
variations  in  the  gravity  field. 

The  nominal  orbital  elements  for  STARLETTE  are  given  in  Table  3. 


Table  3.  Nominal  Orbital  Elements  for  STARLETTE 


Element 

Nominal  Value 

Semi-major  Axis 

Eccentricity 

Inclination 

7, 190  km 

0.0206 

49.83° 

D.  STELLA 

The  STELLA  satellite  was  launched  in  1993  by  CNES,  and  is  virtually  identical  to  the  STARLETTE  satellite. 
As  with  the  LAGEOS  satellites,  both  STELLA  and  STARLETTE  are  spherically-shaped  spacecraft  with 
low  ballistic  coefficients  to  minimize  orbital  uncertainty  caused  by  drag  and  solar  radiation  pressure  forces. 
The  nominal  orbital  elements  for  the  STELLA  satellite  are  given  in  Table  4. 
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Table  4.  Nominal  Orbital  Elements  for  STELLA 


Element 

Nominal  Value 

Semi-major  Axis 
Eccentricity 

Inclination 

7, 178  km 

0.0206 

98.6° 

IV.  Testing  Results 

The  International  Laser  Ranging  Service  (ILRS)  provides  global  satellite  laser  ranging  data  in  order  to 
support  geodetic  research  activities.  Reference  13  contains  further  information  in  regards  to  the  services 
and  capabilities  provided.  The  ILRS  catalogs  SLR  data  to  a  number  of  geodetic  satellites.  These  satellites 
are  typically  design  to  facilitate  the  study  of  Earth’s  gravity  and  contain  laser  retro-reflectors.  This  allows 
highly  precise  laser  ranging  measurements  to  made  which  are  well  suited  for  performing  precision  orbit 
determination. 

SLR  data  from  the  LAGEOS-1,  LAGEOS-2,  STARLETTE  and  STELLA  satellites  will  be  used  to  evaluate 
the  impact  of  using  both  the  state  transition  matrix  calculated  using  only  the  J2  component  of  Earth’s  gravity 
as  well  as  the  state  transition  matrix  calculated  using  the  EGM2008  geopotential  model  to  an  arbitrary  order. 
The  RMS  of  the  residuals  as  well  as  the  number  of  iterations  needed  for  convergence  are  shown.  Additionally, 
differences  in  the  calculated  orbit  ephemerides  are  presented. 

A.  LAGEOS-1  Results 

Table  5  shows  the  RMS  of  the  error  residuals  and  number  of  iterations  needed  for  convergence  for  the  30  day 
WLS-OD  solution  for  the  LAGEOS-1  satellite.  As  can  be  seen,  the  RMS  of  the  error  residuals  are  nearly 
identical  when  either  the  J2  gravity  partial  derivatives  or  the  full  order  gravity  partial  derivatives  are  used. 
However,  the  number  of  iterations  needed  for  convergence  is  reduced  by  two. 

Figure  1  shows  the  RSS  of  the  orbit  position  difference  between  the  two  solutions.  While  the  orbit 
solutions  are  not  identical,  they  only  vary  on  the  order  of  one  millimeter.  This  difference  is  negligible. 

Table  5.  Results  from  LAGEOS-1  30-day  Orbit  Solution 


Method 

Error  Residual  RMS 

Iterations 

J2 

50.95cm 

10 

Full  Order 

50.96cm 

8 

B.  LAGEOS-2  Results 

Table  6  shows  the  RMS  of  the  error  residuals  and  number  of  iterations  needed  for  convergence  for  the  30  day 
WLS-OD  solution  for  the  LAGEOS-2  satellite.  As  can  be  seen,  the  RMS  of  the  error  residuals  are  similar 
when  either  the  J2  gravity  partial  derivatives  or  the  full  order  gravity  partial  derivatives  are  used.  The  orbit 
solution  when  using  the  full  order  gravity  partial  derivatives  does  yield  slightly  lower  error  residuals.  Similar 
to  the  results  from  the  LAGEOS-1  orbit  solution,  the  number  of  iterations  needed  for  convergence  is  reduced 
by  two  when  using  the  full  order  gravity  partial  derivatives. 

Figure  2  shows  the  RSS  of  the  orbit  position  difference  between  the  two  solutions.  Here,  the  orbit 
solutions  vary  on  the  order  of  several  centimeters.  This  difference  is  manifested  primarily  as  an  in-track 
orbit  difference  between  the  two  satellites. 


Table  6.  Results  from  LAGEOS-2  30-day  Orbit  Solution 


Method 

Error  Residual  RMS 

Iterations 

J2 

57.05cm 

11 

Full  Order 

56.98cm 

9 
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C.  STELLA  Results 


Table  7  shows  the  RMS  of  the  error  residuals  and  number  of  iterations  needed  for  convergence  for  the  30 
day  WLS-OD  solution  for  the  STELLA  satellite.  As  can  be  seen,  the  RMS  of  the  error  residuals  is  higher  by 
nearly  25  centimeters  when  using  J2  gravity  partial  derivatives  as  opposed  to  the  full  order  gravity  partial 
derivatives.  Interestingly,  the  number  of  iterations  needed  for  convergence  is  the  same  in  both  cases. 

The  RMS  of  the  error  residual  is  much  higher  here  than  in  either  LAGEOS-1  or  LAGEOS-2  orbit 
solutions.  This  is  due  to  the  nature  of  the  STELLA  orbit,  which  was  designed  to  be  sensitive  to  tidal  effects 
and  temporal  variations  of  the  geopotential.  These  are  often  more  difficult  to  model  over  longer  time  spans. 

Figure  3  shows  the  RSS  of  the  orbit  position  difference  between  the  two  solutions.  Here,  the  difference 
in  the  orbit  solutions  is  much  more  pronounced.  The  RSS  of  the  position  difference  can  vary  by  as  much  as 
3.5  meters.  This  difference  is  manifested  mainly  as  a  cross-track  orbit  difference  (which  accounts  for  about 
2  meters  average  difference),  and  to  a  lesser  degree  as  an  in-track  orbit  difference  (which  accounts  for  about 
1  meter  average  difference). 


Table  7.  Results  from  STELLA  30-day  Orbit  Solution 


Method 

Error  Residual  RMS 

Iterations 

J2 

6.401m 

8 

Full  Order 

6.149m 

8 

D.  STARLETTE  Results 

Table  8  shows  the  RMS  of  the  error  residuals  and  number  of  iterations  needed  for  convergence  for  the  30  day 
WLS-OD  solution  for  the  STARLETTE  satellite.  As  can  be  seen,  the  RMS  of  the  error  residuals  is  higher 
by  nearly  14  centimeters  when  using  </2  gravity  partial  derivatives  as  opposed  to  the  full  order  gravity  partial 
derivatives.  Here,  the  reduction  in  the  number  of  iterations  needed  for  convergence  is  the  most  pronounced, 
requiring  half  the  number  of  iterations  when  using  the  full  order  gravity  partial  derivatives  as  opposed  to 
the  J2  gravity  partial  derivatives. 

Again,  the  RMS  of  the  error  residuals  are  higher  than  either  the  LAGEOS-1  or  LAGEOS-2  orbit  solutions. 
Similarly  to  STELLA,  the  STARLETTE  orbit  was  designed  to  be  sensitive  to  tidal  effects  and  temporal 
variations  of  the  geopotential.  These  are  more  difficult  to  model  over  longer  time  spans. 

Figure  1  shows  the  RSS  of  the  orbit  position  difference  between  the  two  solutions.  Here,  there  is  an 
average  orbit  difference  of  about  1  meter.  This  is  primarily  manifested  as  an  in-track  orbit  difference. 

Table  8.  Results  from  STARLETTE  30-day  Orbit  Solution 


Method 

Error  Residual  RMS 

Iterations 

J2 

4.207m 

12 

Full  Order 

4.071m 

6 

V.  Conclusions 

The  result  show  that  using  the  full  order  geopotential  when  calculating  the  components  of  the  gravity 
gradient  matrix  can  significantly  reduce  the  number  of  iterations  needed  for  the  WLS  orbit  determination 
process  (by  as  much  as  a  factor  of  two).  While  calculating  the  full  order  Hessian  is  more  computationally 
intensive  than  calculating  the  J2  gravity  terms,  the  reduction  in  the  number  of  iterations  still  can  decrease 
the  total  execution  time.  This  is  due  to  the  fact  that  execution  time  for  the  WLS-OD  is  almost  linearly 
related  to  the  number  of  iteration  required. 

Interestingly,  by  using  the  full  order  geopotential  Hessian,  the  RMS  of  residual  error  can  vary  appreciably. 
This  variation  is  much  more  pronounced  for  the  STELLA  and  STARLETTE  satellites,  which  are  in  orbits 
that  are  sensitive  to  higher  order  geopotential  terms.  In  those  cases,  the  solution  orbits  can  differ  significantly. 
Overall,  these  higher  order  gravitational  terms  should  be  used  in  a  WLS-OD  process. 
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Ultimately,  it  is  shown  that  by  using  the  full  order  geopotential  a  substantial  decrease  in  number  of 
iterations  required  for  convergence  can  be  realized  and  a  potentially  better  orbit  solution  may  be  reached. 

Future  work  includes  examining  how  many  higher  order  gravitational  terms  are  necessary  to  realize  these 
gains  and  to  better  characterize  this  method’s  performance  for  satellites  in  different  orbit  regimes.  More 
work  is  needed  to  characterize  the  convergence  behavior  of  both  methods.  It  is  possible  that  by  incorporating 
higher  order  terms  in  the  Hessian  of  the  geopotential  that  radius  of  convergence  in  the  WLS-OD  method  is 
increased. 
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RSS  Position  Difference  (m) 


Figure  1.  LAGEOS-1  RSS  Orbit  Difference  Between  Full  Order  Partials  Solution  and  J2  Partials  Solution 
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RSS  Position  Difference  (m) 


Figure  2.  LAGEOS-2  RSS  Orbit  Difference  Between  Full  Order  Partials  Solution  and  J2  Partials  Solution 
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RSS  Position  Difference  (m) 


Figure  3.  STELLA  RSS  Orbit  Difference  Between  Full  Order  Partials  Solution  and  J2  Partials  Solution 
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RSS  Position  Difference  (m) 


Figure  4.  STARLETTE  RSS  Orbit  Difference  Between  Full  Order  Partials  Solution  and  J2  Partials  Solution 
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